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Bottom baryon masses are calculated based on a 2+1 flavor dynamical lattice QCD simulation. 
The gauge field configurations were computed by the CP-PACS and JLQCD collaborations using an 
improved clover action. The bottom quark is described using lattice NRQCD. Results are presented 
for single and double-b baryons at one lattice spacing. Comparison with experimental values is 
discussed. 



I. INTRODUCTION 

There have been a number of developments since 
our previous systematic study pj of heavy baryons in 
quenched lattice QCD which suggest to us to revisit this 
problem. On the experimental side, masses of four bot- 
tom baryons have been measured since our earlier work. 
On the theory side improved analysis methods should al- 
low for a more precise determination of the masses and 
subsequently a more stringent test of the calculation. As 
well, full dynamical simulations are becoming the norm. 
Dynamical simulations for heavy baryons using lattice 
formulations different from that used in this work have 
been reported in Refs. Q and [jj. 

Besides lattice QCD, heavy baryons have been studied 
in many other approaches. A few recent papers (from 
which the extensive earlier literature may be traced) in- 
clude work on QCD sum rules [3], the quark model 0. 
and the combined heavy quark and l/N c expansions 

With a single b-quark and different flavor and spin 
combinations of up, down and strange quarks eight dif- 
ferent baryons can be constructed. The properties of the 
baryons are summarized in Table [H For the purposes 
of this work exact isospin symmetry is assumed; u and 
d quarks are taken to be degenerate. The experimental 
values for masses of £&, ££, fih and have only become 
available in the past yearfj, flQ. [Tl|. 

A variety of approaches are being pursued to incor- 
porate dynamical quark effects in lattice QCD simula- 
tions. Ideally one would like to have u/d and strange 
(2+1 flavor) dynamical effects. As well, the light quarks 
(u/d) should have small masses and the lattice volume 
should be large. It is clear that these requirements stress 
the computing resources available to even the largest lat- 
tice QCD collaborations so some compromises have to be 
made. It is not feasible for us to generate our own dynam- 
ical gauge field configurations. Instead 2+1 flavor dy- 
namical configurations, made available by the CP-PACS 
and JLQCD collaborations[I2,[l3j], were used. These are 
based on the Iwasaki RG gauge field action[3| and an 
improved clover actionpjl for the quarks. 

Lattice NRQCD[l(| is used for the heavy quark. For 
the lattice spacing considered here, only the b quark is 
heavy enough for its mass to lie above the cut-off scale. 



Table I: Properties of single-b baryons showing valence con- 
tent (q=u/d), spin parity, isospin and mass (in GeV). The 
quantity s; is the total spin of the light quark pair. 
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Table II: Properties of double-b baryons showing valence con- 
tent (q=u/d), spin parity and isospin. 



Baryon 


quark content 


J p I 


E*6 
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1+ 1 
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Charm is too light to be simulated by NRQCD and so, 
unlike Ref. charmed baryons are not considered in 
this study. 

The parameters in the calculation, namely u/d, s and 
b quark masses as well as the overall scale, are fixed by 
calculations done in the meson sector. Masses of heavy 
baryons are then predictions of the simulation. In addi- 
tion to single-b baryons listed in Table U we also calculate 
the masses of the double-b baryons with properties given 
in Table HH At present there are no data for double- 
b baryons but it is hoped that eventually they will be 
observed in some future experiment. 

In Sec.|n]details of the simulation are presented. There 
is a brief summary of the actions and the lattice param- 
eters. The analysis method is also discussed. 
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Table III: Lattice parameters. The lattice size is 20 3 x 40 
with Iwasaki gauge action /3 = 1.90 and quark action clover 
coefficient csw = 1.715. The masses of u/d and s quarks are 
encoded in the hopping parameters K q and n a respectively. 



configurations 




Kg 


U 


480 


0.1358 


0.1358 


0.8396 


576 


0.1364 


0.1358 


0.8415 


576 


0.1368 


0.1358 


0.8425 


576 


0.1370 


0.1358 


0.8432 


480 


0.1358 


0.1364 


0.8405 


576 


0.1364 


0.1364 


0.8422 


576 


0.1368 


0.1364 


0.8433 


576 


0.1370 


0.1364 


0.8439 



Results are given in Sec. Hill and compared to available 
experimental data. Limitations of this work and future 
directions are discussed in Sec. IIVI 



II. NUMERICAL SIMULATION 
A. Simulation details 

Gauge field configurations incorporating the full dy- 
namical vacuum polarization effects of u,d and s quarks 
are employed in this work. These gauge configura- 
tions were generated by the CP-PACS and JLQCD 
collaborations [13, [13] and made available through the 
Japan Lattice Data Gridfl7l]. 

The gau ge fi eld action is an improved action developed 
by Iwasakifl4] and in addition to the standard plaquette 
term it contains six-link operators with coefficients tuned 
so that the action lies close to the renormalized trajec- 
tory. The quark action is of the clover type[15|. The 
coefficient of the clover was tuned nonperturbativelvfl8| 
and is the same for all values of sea-quark mass. Hy- 
brid Monte Carlo was used for the two light flavors (u/d) 
and the Polynomial Hybrid Monte Carlo was used for the 
strange sea quark[19]. 

In this work, a subset of the CP-PACS/ JLQCD 2+1 
flavor configurations at /3 = 1.90 with lattice size 20 3 x 40 
is used. The coefficient of the clover term in the quark 
action was 1.715 (see Ref. [jji])- Other lattice parameters 
are given in Table HTT1 The quantity Uq equals 1/3 of the 
trace of the mean gauge-field link in Landau gauge and 
is used as the tadpole factor in the NRQCD action. 

A non-relativistic action[l(| is used to describe the b 
quark. This approach is useful when the bare mass of the 
quark is larger than the cutoff scale. The particular form 
of the action used in this work is essentially the same 
as that used in our previous works [H Terms up to 
0(1/ Mq) in the heavy quark mass M are retained. The 
only difference is that while [l|, [2(J used an anisotropic 



lattice, here the lattice is isotropic. The details of the 
action are given in the Appendix. 

Simulations were done for three values of the heavy- 
quark bare mass: 2.28, 2.34 and 2.40 in lattice units. This 
allowed for interpolation to the physical b-mass value. 

Correlation functions were calculated using local 
hadron operators at the source and sink. For heavy 
baryons the operators used were exactly the same as 
detailed in Ref. [l|. In addition to baryon correlators, 
correlators were calculated for pseudoscalar and vector 
mesons in the light sector (u/d,s), the heavy-light sector 
(B-mesons) and the heavy-heavy sector (T). The meson 
calculations were used to determine the "physical" val- 
ues for the quark masses and the overall scale. As well 
they provide some additional predictions which test the 
calculation. 

Non-relativistic quarks propagate only forward in time. 
For combining relativistic quarks with non-relativistic 
quarks it is convenient to use Dirichlet boundary condi- 
tions for the light quark propagators; propagation across 
the time boundary is not allowed. If the source time is 
set in from the time boundary (four time steps is used in 
this calculation) it can be verified that meson masses in 
the light sector are the same, within statistical errors, as 
those computed with more usual periodic time boundary 
conditions. 

In order to get some extra suppression of statistical 
fluctuations and reach the level of precision that we would 
like, multiple sources were used. For each gauge field con- 
figuration a set of correlators was calculated using differ- 
ent space-time points as the source point. This set was 
averaged and the average value was used as the represen- 
tative correlation function for the configuration. Eight 
sources per configuration were used in this calculation. 



B. Analysis 

Correlation functions were fit with a sum of exponen- 
tials 

n-l 

i=0 

over a fixed time range which for our standard analysis 
extended for 27 time steps starting one time step past 
the source position. For mesons three exponential terms 
were found to be adequate while for baryons four terms 
were used. The advantage of a multi-exponential fit over 
methods such as fitting over an "effective mass" plateau is 
that it is less subjective. Also it makes better use of the 
correlation function data at times where the statistical 
error is small. 

To stabilize the fits a constrained fitting method|21] 
was used. The usual x 2 is augmented with a term which 
acts to prevent the fit parameters from straying outside 
some sensible (but broad) range. The constraint term for 
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the coefficients Zi was taken to be 

gfc^. (2) 

The priors were chosen to be Zi = o gi = g(0)/n, which 
is a very loose constraint. The constraint term for the 
exponents (energies) has a similar form 

E^y^- (3) 

i=0 Ei 

One might be tempted to choose a constraint which is 
minimally biased such as equal spacing of the energies 
Ei — Ei-i = SE and erg. — SE. In practice it was found 
that to prevent the occurrence of obviously spurious solu- 
tions, where for example two terms have the same expo- 
nent and coefficients of opposite sign, a somewhat tighter 
constraint is needed. For our fits we use E\ — Eq = 0.5 
and Ei—Ei-i = 1.0, i > 1 with Eq equal to 0.5 for mesons 
and 0.8 for baryons. The were chosen to be propor- 
tional to the spacing between the Ei with a common re- 
duction factor taken to be 0.83. With this setup spurious 
solutions were avoided and the ground state energy was 
very stable with respect to changes in the constraint pa- 
rameters. For example, it was found that changing prior 
parameters by 20% led to changes in single-b baryon sim- 
ulation energies (which are most sensitive) of about 1%. 
This effect is incorporated into the estimated systematic 
uncertainty. 

To determine the statistical errors, bootstrap analy- 
sis was used. A sample of 600 bootstrap ensembles was 
created and a complete analysis was carried out for each 
ensemble. In this way the uncertainty in determining the 
quark mass parameters and the scale is incorporated into 
the bootstrap error estimate of the final result. 

A few systematic effects were examined explicitly. 
These include sensitivity to the choice of time range for 
fitting the correlators, to the determination of the fa- 
quark mass and to the choice of extrapolation function 
in u/d and s quark mass. In addition, there are system- 
atic errors associated with the omitted higher order ef- 
fects in the NRQCD action. These are discussed in great 
detail by Gray et al. [HI] for the T system. For heavy- 
light hadrons the appropriate power counting is different 
than for quarkonium[23j . The expansion parameter is 
p/Mq where the typical momentum of the heavy quark 
p ~ Aqcd, independent of quark mass. Radiative cor- 
rections to the action then are 0{a B KQCD / Mq) relative 
to the leading kinetic term[23] which we estimate can 
induce a 1.5% relative uncertainty in the simulation en- 
ergies. The NRQCD action is corrected for 0(a 2 ) lattice 
spacing errors at tree level. The leading discretization 
corrections are 0(a s (aAQCD) 2 ) and a 0.5% fractional 
systematic error is assigned to these effects. The action 
used in this work includes 0(1/ Mq) terms and higher 
order relativistic effects are assumed to be negligible. 



III. RESULTS 
A. Mesons 

The first task is to determine the quark mass parame- 
ters and overall scale. We start with the u/d and s sec- 
tor and use the pion, rho meson and phi meson masses 
as experimental input. Lattice masses for these mesons 
are calculated with the eight ensembles in Table III. For 
each correlator the valence mass was taken to be equal 
to the u/d mass or to the strange mass. No partially 
quenched correlators, with a valence mass different from 
a sea-quark mass of the same flavor were used in this 
work. 

The lattice meson masses were then fit as a function of 
quark mass using the vector Ward identity (VWI) defi- 
nition of the quark mass m = (1/k — l//t cr )/2 where K cr 
is the point where the pseudoscalar meson mass vanishes 
when all valence and sea quark masses have this hopping 
parameter. 

The fitting functions are motivated by the study of 
light quark masses in Ref. [13]. For the pseudoscalar 
meson (with degenerate valence quarks) it was found that 
the three terms 

bim v + 62(2m 9 + m s ) + bzm v (2m q + m s ), (4) 

where m v , m q and m s are valence, u/d and s quark 
masses from the VWI definition, gives a good descrip- 
tion of the mass squared lattice data. Note that 
free parameter in the fit as are the coefficients 6i, 62, 63. 
For the vector meson the mass was fit using 

ci + c 2 m v + c 3 (2m q + m s ) + c±ml (5) 

which takes into account a slight nonlinear dependence 
on the quark mass. In determining the coefficients of this 
fit fixed to the value obtained in the pseudoscalar 

meson fit. 

After the fit parameters are determined one can solve 
for the values of the hopping parameters and the lattice 
spacing that reproduce the input experimental numbers. 
The results are n q = 0.13784(2) and n s = 0.13618(9) with 
a value of the inverse lattice spacing a -1 = 1.89(3)GeV. 
We note that the inverse lattice spacing obtained here is 
slightly different from that quoted in [13j but it has to 
be remembered that the number of configurations, the 
number of ensembles and the details of the mass extrap- 
olation are also different. 

The b-quark mass is determined using the T mass as 
input. In NRQCD the quark mass has been removed 
from the action. The zero-momentum hadron correla- 
tor yields a simulation energy to which the renormalized 
quark mass and energy shift must be added to get the 
hadron mass. Alternatively the hadron mass can be de- 
termined from the kinetic energy. For this purpose the 
correlators for bb mesons carrying a unit of momentum 
were calculated. The meson mass M (lattice units) is 
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obtained using 



M = 



2n 2 



N*(Es im (p) -E aim (0)) 



(6) 



where N s is the spatial extent of the lattice, p — 2ir/N s 
and the difference in simulation energies is just the kinetic 
energy. 

Kinetic energies were calculated for T at all the heavy 
bare quark masses listed in Sec. Ill Al and for all eight 
ensembles. It was found that kinetic energies were prac- 
tically independent of the values of the sea quark masses 
with no systematic trend which would allow for extrap- 
olation. Therefore for each ensemble an independent de- 
termination of the bare b-quark mass was made by in- 
putting the experimental T mass and the inverse lattice 
spacing obtained above. The values (in lattice units) 
varied in the small range from 2.375(78) to 2.416(77). 
A value in the middle of this range 2.391(78), obtained 
from the ensemble K q = 0.1368, k s = 0.1364, was used 
as our nominal value for the b mass. The maximum and 
minimum values were used to estimate a systematic un- 
certainty in the hadron masses. The effect of choosing 
different values for the b-quark mass was found to be 
small, typically changing the final hadron mass by an 
MeV or so. 

Having determined the physical point for the b-quark 
mass one can get the splitting between the T and the iji, 
just from interpolating the difference in simulation en- 
ergies for the 66 vector and pseudoscalar channels. The 
calculated T — rjb mass difference is given in Table IV and 
is smaller than the recent value reported by the BaBar 
collaborationjH]. The underestimate of the quarkonium 
spin splitting has been a common feature of simulations 
done with 

NRQCDJUStIHI]- The systematic study by 
Gray et al. [22] (see their Fig. 14) shows that this quantity 
is quite sensitive to the continuum extrapolation for light 
sea-quark masses. Gray et aL[22] predicted a continuum 
extrapolated value of 61(14)MeV, somewhat larger than 
our value at a non-zero lattice spacing. Note that part 
of the difference between their value and ours is due to 
0(v 6 ) terms in the NRQCD action which we include and 
they do not. From a test run on a single ensemble it is 
estimated that these terms decrease the spin splitting by 
about 15%. 

Although the main motivation for this work was to 
study heavy baryons, the masses for heavy-light (B) 
mesons were also calculated. The masses were com- 
puted relative to the T which was an input to the cal- 
culation. Above, the kinetic energy was used to deter- 
mine the T mass but alternatively it is given as M*r = 
Eg iT n + 2(ZM — E shift) where Z is the mass renormal- 
ization factor and E s hift 15 an additive mass shift. The 
quantities Z and E s hi ft are independent of the hadronic 
state[24| so for a hadron, such as the B-meson, contain- 

ing a single b quark Mb = E y sin [ + ZM - E shi f t . The 
mass can then be obtained using 



Table IV: Lattice results for meson masses and spin splittings 
(in GeV) compared to experimental values. The T — rjb result 
is from [25l |. other experimental values are from Q. 



mass-M T /2 



Experiment 



T — rjb 
B 
B" 
B s 
Bt 



0.527(6) (8) 
0.530(7)(J°) 
0.617(3)(10) 
0.658(4)(11) 



B* 
Bt - 



- B 
B s 



0.039(1)(7) 
5.257(6)(8) 
5.300(7) (J ) 
5.346(3)(10) 
5.388(4)(11) 
0.043(3)(|) 
0.042(2)(j) 



0.07l4(i)(27) 
5.2793(4) 
5.325(1) 
5.366(1) 
5.412(1) 
0.0458(4) 
0.0461(15) 



The results, extrapolated in light quark mass and con- 
verted to physical units, are shown in Table IIVI along 
with the experimental values from the PDG[8]. The first 
error in the lattice results is the statistical (bootstrap) 
error. The second error incorporates the changes that 
result from changes in the time range used in fitting the 
correlators, the uncertainty in determining the b-quark 
mass and the choice of u/d and s quark mass fitting func- 
tion. For the mass range used in this present simulation 
the heavy-light meson masses are consistent with a linear 
dependence on quark mass 



Ci + c 2 m v + c 3 (2m q + : 



(8) 



J sim 



1 



(M- 



E m ). 

sim J 



(7) 



where m v ,m q and m s are the VWI quark masses. This 
was the form used to extrapolate for u/d and interpolate 
for s. The systematic error incorporates an estimate of 
sensitivity to non-linear quark mass dependence made by 
adding a quadratic term as in Eq. |(5]). Also included in 
the systematic error are estimated uncertainties associ- 
ated with radiative and discretization corrections to the 
NRQCD action. 

The heavy-light meson results look quite reasonable al- 
though there is a systematic underestimate of the mass 
by about 25MeV compared to experimental values. It is 
plausible that a common effect underlies this trend. At 
this stage, we don't know the effect of changing the lat- 
tice spacing so it is natural to suspect that the overall dis- 
crepancy in Table |IV] is due to a small uncorrected lattice 
spacing error. Assuming the lattice spacing error associ- 
ated with the light quark action is C(o 2 Aq CD ) one gets 
the estimate a 2 AQ C . £) (mass-Mx/2) ~ 0.015GeV which is 
roughly the size of the observed discrepancy. Note that 
for heavy-light mesons, where the wavefunctions are not 
so strongly affected by short distance interactions as in 
bottomonium, there is good agreement of the calculated 
spin splittings with experimental values. 



B. Baryons 

The calculation of baryon masses proceeds very much 
like that for heavy- light mesons. The simulation ener- 
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Figure 1: A Ej correlator as a function of lattice time. The 
lattice boundaries are at t equals 1 and 40. The source is at t 
equals 5. The fit (solid line) is done including points t equals 
6 to 32. 



Table V: Lattice results for masses of single-b baryons (in 
GeV) compared to experimental values. 



mass- M T /2 



Experiment Ref. 



A 6 
St 



0.911(21)(jl) 
1.065(16)(i) 
1.112(26)(3) 

1.276(10)(8) 
1.314(18)(1?) 

1.051(17) (S) 
1.173(12X8) 
1.220(21)Gj) 



5.641(21)$) 
5.795(16) (£) 
5.842(26)(fg) 
6.006(10)(iB) 
6.044(18)(1?) 
5.781(17) (g) 
5.903(12)(i|) 
5.950(21)^?) 



5.620(2) 
5.8115(30) 
5.8327(34) 
6.165(16) 



[8] 
[9] 
[9] 
[10] 



5.7929(30) [11] 



gies were determined by doing four-term exponential fits 
to the baryon correlation function. A typical example, 
showing the quality of the simulation data and of the fit, 
is given in Fig. [TJ For single-b baryons, masses were cal- 
culated using the baryon analog of Eq. |(7|). The resulting 
masses were first interpolated in the b-quark mass to the 
physical point determined by My. Then the u/d and s 
quark mass dependence was fit using Eq. J8]). A linear 
quark mass dependence was consistent with all simula- 
tion data. A quadratic valence mass dependence was 
used to estimate a systematic uncertainty. 

The results for single-b baryons are tabulated in Ta- 
ble [V] and plotted in Fig. [2] with statistical errors and in 
Fig. [3] with combined statistical and systematic errors. 
The As and E& show significant sensitivity to the inclu- 
sion of a quadratic term in the mass fit; other masses 
are not changed very much by this term. The present 
results show a vast improvement in statistical precision 
compared to our previous study [lj]. Multi-exponential 
constrained fitting which uses time-correlation informa- 
tion over a large range, including times where the corre- 



5.6 



A, 



a, 



"b \ 



Figure 2: Masses of single-b baryons. The diagonally-hatched 
boxes are lattice results with statistical errors only. Solid bars 
(red) are experimental values. 




Figure 3: Masses of single-b baryons. The diagonally-hatched 
boxes are lattice results with combined statistical and system- 
atic errors. Solid bars (red) are experimental values. 



lation function statistical errors are small, plays an im- 
portant role in this improvement. Also shown in Table 
IVl is the quantity mass-Mr /2. This is the actual quan- 
tity that the lattice simulation provides and, since most 
of the heavy hadron mass is due to the b-quark mass, 
it is a rough measure of the quark and gluon interaction 
energy. With the present analysis methods this quantity 
is determined to a few percent. 

With the exception of fib the results are in good agree- 
ment with the experimental values. However, one should 
not overinterpret this agreement since scaling (with lat- 
tice spacing) has not been checked and effects of using 
more realistic u/d quark masses still have to be consid- 
ered. 

The large discrepancy between our calculated fib mass 
and the value reported in [l(J is perplexing. To under- 
stand how puzzling it is, consider the basic physical idea 
behind heavy quark effective theory: a heavy quark acts 
essentially as a static color source and so mass differences 
between states with different light-quark configurations 
should be independent of heavy quark mass (up to cor- 
rections inversely proportional to the heavy mass). This 
suggests a direct data-to-data comparison of mass dif- 
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Table VI: Lattice results for masses of double-b baryons (in 
GeV). 





mass- My 


mass 


Hfcb 
i — * 
— bb 

Qbb 

^66 
—66 — —66 
Qlb ~ ^66 


0.667(13)$) 
0.691(14)®) 
0.762(9)(}|) 
0.786(10)61) 


10.127(13) (i) 
10.151(14)®) 
10.225(9)(}|) 
10.246(10)(i|) 
0.026(8)®) 
0.025(7) (J 1 ) 



Figure 4: Experimentally measured values of masses of single- 
charm (dashed, blue) and single-bottom baryons (solid, red) 
relative to the lowest lying state A. 
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Figure 5: Masses of single-b baryons from Jenkins[7] 
(vertically-hatched boxes, blue). The diagonally-hatched 
boxes are lattice results and solid bars (red) are experimental 
values. 



ferences in single-charm and single-bottom baryons as 
shown in Fig. [H The measured masses of the singly- 
heavy S, E* and 3 baryons fit the expected pattern but Q 
shows a large discrepancy, the same behavior as observed 
with our lattice simulation. The application of heavy 
quark effective theory to singly-heavy baryons was for- 
malized by Jenkins [29]. In this work a combined expan- 
sion in the inverse of the heavy-quark mass, in l/N c and 
in SU(3) flavor symmetry breaking was carried out. Mass 
formulas were derived which then allow some masses to 
be predicted in terms of other experimentally measured 
masses. This is a more rigorous version of our data-to- 
data comparison. The updated predictions for single-b 
baryons from JenkinsQ are shown in Fig. [5] along with 
our lattice results and the experimental values. There 
is consistency between our lattice results and the effec- 
tive theory analysis of 0, H§] across the whole spectrum 
which further accentuates the fib puzzle. 

The masses of double-b baryons are calculated using 



M bb - E^t, 



M T - E 



(T) 



(9) 




Figure 6: Lattice results for masses of double-b baryons show- 
ing statistical errors only. 



The final results are listed in Table I VII and shown in 
Figs. [6] and [7] with statistical and combined errors re- 
spectively. As yet there are no experimental values to 
compare with these calculations. The mass difference be- 
tween spin 3/2 and spin 1/2 states in the doubly heavy 
sector is interesting because in the quark model and 
in heavy quark effective theory the baryon spin split- 
ting can be related to the spin splitting for heavy-light 




The masses are interpolated and extrapolated as above. 



Figure 7: Lattice results for masses of double-b baryons show- 
ing combined statistical and systematic errors. 
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mesons [33, HH, [H, Simple arguments suggest that 

3 

AMharyon ~ -^AM meson . (10) 

Within errors, our lattice spin splittings are consistent 
with this relation but are not precise enough to test it 
very stringently. 



sea and valence quarks are treated using different actions 
needs further consideration. 

In our previous papers charmed baryons were also 
studied[l|, [52] and one would like to test one's capability 
of doing calculations in this sector. At the lattice spacing 
used here the charm quark is too light for NRQCD so the 
treatment of charmed baryons is left for future work. 



IV. SUMMARY AND DISCUSSION 



APPENDIX 



In this paper the masses of bottom baryons were calcu- 
lated in a 2+1 flavor dynamical lattice simulation. The 
clover action was used for light (u/d and s) quarks and 
the b quark was described by NRQCD. Using constrained 
multi-exponential fitting for hadron correlators led to 
much more precise results than our previous studies. Sin- 
gle bottom baryons whose masses can be compared to 
experimental values show good agreement with the ex- 
ception of the recently measured fib. 

The fib is puzzling since our lattice result is consistent 
with ideas based on heavy quark effective theory which 
allow for a prediction of the fib mass using only empirical 
input. 

Predictions are made for the still unobserved single- 
b baryons and for the double-b baryons. Experimental 
observation of any of these states would be extremely 
interesting and might shed some light on the fib puzzle. 

This study should be considered as the first step in a 
program to do precision calculations in the heavy baryon 
sector. To go further, some systematic improvements 
have to be made. Radiative corrections to the NRQCD 
action are one of the major contributions to the system- 
atic error and have to be dealt with. The light quark 
action is not corrected for lattice spacing errors to the 
same extent as the heavy quark action. Calculations 
are needed at more than one lattice spacing to enable 
a continuum extrapolation and to estimate reliably lat- 
tice spacing errors. 

The present simulation is done in a region where the 
u/d quark masses are about 0.4 times the strange mass 
and some baryon masses exhibit a sensitivity to light 
quark mass extrapolation. Smaller values for the u/d 
mass would be highly desirable to insure that the simu- 
lation captures more accurately the dynamical sea quark 
effects and that the extrapolation to the physical point 
can be put on a firm theoretical basis. Using an array 
of algorithmic improvements, the PACS-CS collaboration 
has produced configurations with the clover action push- 
ing the light quark masses to near the physical region (34|. 
As of this writing, these configurations are not avail- 
able for our use. Other fermion formulations could be 
considered, such as staggered fermions[35|], domain-wall 
fermions[3a| or twisted mass Q CD (3?J which presently 
operate at u/d to s-quark mass ratios as small as 0.1, 
0.217 and 0.16 respectively. However, how to combine 
NRQCD with such approaches may require some atten- 
tion. Also, the possibility of a hybrid calculation where 



The heavy quark action is described using NRQCD [Hi . 
The heavy quark propagator is given by 



Gr+i = 1- 



clHb \ ( ', aHA\ n U\ 



xl- 



2 

aH A 
2n 



1 - 



1 - 



2n J U 
aH B 



G T , (Al) 



with n = 5 used in this work. The Hamiltonian is sepa- 
rated into two terms, H — Ha + Hb, with Ha containing 
the kinetic piece Hq and the term proportional to C10 (de- 
fined below). 

The Hamiltonian contains all terms up to 0(1/Mq) in 
the classical continuum limit: 



H = H + SH } 
-A( 2 ) 

Hn = 



2Mn 



(A2) 
(A3) 



SH = <5iJ (1) + 8H^ +O(l/A/ 4 ), (A4) 



U£ 2 AT 



c 5 - 



24M ' 



(A5) 



SH 



( 2 ) _ c 2 ig 



C3 9 



C/ n 4 8M n 2(A ' E " E ' A) r^^r2 



C/ 4 8M 2 



xer • (A x E - E x A) - c 6 



»(Af 



2)\2 



16nM 2 



■,(A6) 



SH w (aw) 2 r A(2) 

c 9 ig 2 /ExE BxB\ 



8A/ 3 



U§ C/ 8 



c w g 2 ( E 2 B 2 \ a 2 (A( 2 )) 3 



8M 3 \U* U* 



192n 2 M 3 



(A7) 



The tildes indicate that the leading discretization errors 
have been removed. In particular, 



Ei — 

1 

2 



(A8) 
(A9) 



where 
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+Ul(x - fi)F^(x - fi)U^{x - /i) 

-(M <- v)] + I (± - lj F^x). (A10) 

The last term in F^ v {x) corrects for the fact that the 
gauge field link multiplied by a tadpole factor is no longer 
unit ary [3^ . 

The spatial lattice derivatives are given by 



aAjG(x) 
aA^G(x) 

a 2 A^G(x) 



^L[Ui(x)G(x + l) -U}(x- i)G(x - 1)}, 

(All) 
(A12) 



Ui(x) 



G(x + i)-G(x), 



G(x) - U ^ X TT — G(x - i) 



Uijx) 

, uUx-i) 



U 

G(x + 1) - 2G(x) 
G{x-i), 



U, 



(A13) 



(A14) 



A, 

A (2) 

A(2) 

A (4) 



G 

(2) 



A< 2 »--A( 4 > 
12 

EN 25 " 



(A15) 
(A16) 

(A17) 
(A18) 
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